A threshold induced phase transition in the kinetic exchange models 
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We study an ideal-gas-like model where the particles exchange energy stochastically, through 
energy conserving scattering processes, which take place if and only if at least one of the two 
particles has energy below a certain energy threshold (interactions are initiated by such low energy 
particles). This model has an intriguing phase transition in the sense that there is a critical value 
of the energy threshold, below which the number of particles in the steady state goes to zero, and 
above which the average number of particles in the steady state is non-zero. This phase transition is 
associated with standard features like "critical slowing down" and non-trivial values of some critical 
exponents characterizing the variation of thermodynamic quantities near the threshold energy. The 
features are exhibited not only in the mean field version but also in the lattice versions. 
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I. INTRODUCTION 

The kinetic theory of gases had played a pivotal role in 
the development of statistical mechanics, which is more 
than a century old. This theory describes a gas as a col- 
lection of a large number of particles (atoms or molecules) 
which are constantly in random motion, and these rapidly 
moving particles constantly collide with each other and 
exchange energy. In the ideal gas, this energy is only ki- 
netic. Recently physicists have been studying two-body 
kinetic exchange models in the socio-economic contexts 
in the rapidly growing interdisciplinary field of "Socio- 
physics" [H and "Econophysics" 2]. The two-body ex- 
change dynamics has been developed in the context of 
modeling income, money or wealth distributions in a so- 
ciety I3l4l oil ■ and modeling opinion formation in the so- 
ciety |ll| - [l4| . analogous to the kinetic theory model of 
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ideal gases. These studies have given deeper insights and 
different perspectives in the simple physics of two body 
kinetic exchange dynamics. In this context of wealth ex- 
change processes, Iglesias et al. [l^ [1^1 had considered 
a model for the economy where the poorest in the so- 
ciety (atom with least energy in the gas) at any stage 
takes the initiative to go for a trade (random wealth / 
energy exchange) with anyone else. Interestingly, in the 
steady state, one obtained a self-organized poverty line, 
below which none could be found and above which, a 
standard exponential decay of the distribution (Gibbs) 
was obtained. 



Here, we study a model where TV particles, interact 
among themselves through two-body energy {x) conserv- 
ing stochastic scatterings with at least one of the parti- 
cles having energy below a threshold (poverty line in 
the equivalent economic model) . The states of particles 
are characterized by the energy {xi\, i = 1,2, ...,A^, 
such that Xi > 0, Vi and the total energy E = '^^Xi 
is conserved (= N here, such that the average energy of 
the system E — E/N =1). The evolution of the system 
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is carried out according to the following dynamics: 

x'^ = {l-e){x<+x,), (1) 

where xf < 9 (threshold energy or " poverty line") and 
e (0 < e < 1) is a stochastic variable, changing with 
time (scattering). It can be noticed that, the quantity x 
is conserved during each collision: xf + x'^ — xf + Xj . 
The question of interest is: "What is the steady state 
distribution p{x) of energy x in such systems?" 

In the standard case, when the threshold energy goes to 
infinity {0 — ^ 00), we know that the steady state energy 
distribution will be the exponential Gibbs distribution 
{p{x) ~ exp(— cc)) 1^. However, when a finite threshold 
energy is introduced {0 > 0), several new and intriguing 
features appear. These features are exhibited not only 
in the mean field version (with infinite-range interaction, 
pairs of particles randomly chosen from N particles) but 
also in the lattice versions (with nearest neighbor inter- 
actions, i.e. exchanges between the nearest neighbors on 
lattice sites). 

II. MODEL SIMULATIONS AND RESULTS 

A. The Model 

We simulate a system of N particles (agents). At any 
time we select randomly a particle i. If the energy 
of the particle is below a prescribed threshold energy 9, 
then it collides with any other random particle j (in the 
mean field model) which can have any energy whatsoever, 
and the two particles will exchange energy according to 
the Gibbs-Boltzmann dynamics of Eq. ([T]). After each 
such successful collision, the time is incremented by unity. 
The dynamics will continue for an indefinite period, un- 
less there is no particle left below the threshold energy, 
in which case the dynamics will freeze. If the dynamics 
gets frozen (when Xi > 6 for all i), we employ a 'mild' 
perturbation such that a randomly chosen particle will be 
dropped to the lower level (< 9) by giving up its energy 
to anyone else (to ensure total energy conservation). It 
can be shown that the addition of this perturbation does 
not alter the relevant quantities for a thermodynamically 
large systems, and simply ensures ergodicity in the sys- 
tem. After sufficiently large time t > t, a steady state 
is reached when the energy distribution p{x) (and also 
other average quantities) do not change with time. We 
start with different initial random configurations, where 
the states of particles are characterized by the energies 
{xi}, i = 1,2, . . . , N , which are drawn randomly from 
an uniform distribution such that Xi > 0, \/i and the 
average energy E = J2i is set to unity. We find 

the system to be ergodic (the steady state distribution 
p{x) is independent of the initial conditions {xi}), and 
we take steady state averages over all such independent 
initial conditions to evaluate the quantities of interest. 
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FIG. 1. (Color online) Energy distribution p{x) in the steady 
state (t > r), for different values. Th inset shows semi-log 
plot of the energy distribution. The tail of the distribution 
is Gibbs-like (A'' = 10^; Mean field model with average taken 
over many independent initial conditions). 



We study mainly three cases (a) mean field (or infinite 
range) case where i and j in Eq. ((1} can represent any two 
particles/agents in the system; (b) one dimensional case 
where j = i H along a chain and (c) two dimensional 
case, where j = i ± S where S represents neighbors of i. 
In our studies we consider a 2D-square lattice. 

We observe that for finite values of the energy thresh- 
old 9, the steady state energy distribution is no longer 
the simple Gibbs-Boltzmann distribution. We also find 
that O (= Jq p{x)dx), the average number of particles 
below the threshold energy in the steady state, is zero 
for 9 values below or at a critical threshold energy 9c, 
and for 9 > 9c, O is non-zero. The steady state value of 
O, the average number of particles below the threshold 
energy 9 is seen to act like an "order parameter" of the 
system. We study the relaxation dynamics in the sys- 
tem: the relaxation of 0{t) to the steady state value of 
O (= 0{9) for t > t{9), the "relaxation time"). We find 
t{9) grows as 9 approaches 9c, and eventually diverges 
at 9c- The details of the results are given below. 



B. Results: Mean field model 

In the mean field (long range) model, we first look 
for any particle (z) with energy Xi < 9 and then this 
particle is allowed to interact with any other particle (j), 
following Eq. This continues until either the steady 
state, or a frozen state with Xi > 9 for all i is reached. 
In the case of frozen state, as mentioned earlier, any one 
particle is picked up randomly and it loses its energy to 
any other (randomly chosen) particle, and goes below the 
threshold 9. This induces further dynamics. Eventually 
steady state is reached. We study this steady state energy 
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FIG. 2. (Color online) Simulation results for the variation 
of O, the average number of particles below the threshold 
energy 6 in the steady state {t > r), against threshold energy 
e. (Inset)(a) Shows the results at O for 61 = 0.59(< Oc) 
as TV -> oo; (b) Shows scahng fit {6 - Ocf with P ~ 0.97. 
(A'' = 10""; Mean field model with average taken over many 
independent initial conditions). 



distribution p{x) (see Fig. [T]), and the order parameter 
O = Jg p{x)dx (see Fig. [5]), showing a "phase transition" 
at Oc ~ 0.607±0.001. A power law fit O ~ {O-Ocf gives 
13 ~ 0.97 ±0.01. 

We also studied the relaxation behavior of O. At 9 = 
0c, the 0{t) variation fits well with r*; S ~ 0.93 ± 0.01 
(see Fig. The relaxation time r is estimated numer- 
ically from the time value at which O first touches the 
steady state value 0{9) within a pre-assigned error limit. 
We find diverging growth of relaxation time r near 6 = 6c 
(see Fig. 4), showing 'critical slowing down' at the criti- 
cal value 9c- The values of exponent z for the divergence 
in T \0 — Oc\~^ have been estimated (for both > 9c 
and < 9c)- For the mean-field model, the fitting value 
for exponent z ~ 0.83 ± 0.01. 

We have also studied the universality of this behavior 
by generalizing the dynamics in Eq. ([1]) to 



(l-ei)a;<+(l-e2)a 



(2) 



where ei and 62 are the random stochastic variables 
within the range [0,1]. The critical point 9c shifts to 
9c ~ 0.69 {9c ^ 0.61 for ei= £2 = e ). The transition be- 
haviour is seen to be universal near the critical point 9c, 
but the critical point depends specifically on the model 
(see Fig. [S]). 



C. Results: One dimensional model 

In one dimensional lattice version, the particles are ar- 
ranged on a periodic chain. At any time t, we randomly 
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FIG. 3. (Color online) Variations of O versus time t, shown 
for different values. At critical value 9c, order parameter 
follows a power law decays with exponent 5 ~ 0.93 (shown 
by the solid line). {N = 10^; Mean field model with average 
taken over many independent initial conditions). 
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FIG. 4. (Color online) Variation of r versus 9. (Inset) Scaling 
fit r '--^ 1^ — 9c\^^ , with exponent z ~ 0.83 (Mean field case; 
iV = 10^). 



select a lattice site i. If the energy of the corresponding 
particle is below a prescribed threshold energy 9, then 
it collides with any one randomly chosen nearest neigh- 
bors j{— i ± 1) which can have any energy whatsoever, 
and the two particles will exchange energy according to 
Eq. ([T]). After each such successful collision, the time 
is incremented by unity. This process is continued until 
steady state is reached. The steady state order param- 
eter O variations against theshold 9 is shown in Fig. [SI 
with exponent l3 ~ 0.41 ± 0.02 and 9c ~ 0.810 ± 0.001. 
The fitting value for exponent z turn out to be around 
1.9 ± 0.05 (see Fig. El). Also we find 5 ~ 0.19 ± 0.01 (see 
Fig. [5]) . It may be noted that a recent study of a related 
chain model with such energy cut-off for kinetics, where 
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FIG. 5. (Color online) Variation of steady state order param- 
eter O{0) against 6 for dynamics following Eq. U]) (denoted 
by red squares) and Eq. ((2]) (denoted by green circles). (Inset) 
Shows O vs. [9 - Oc) for both cases. (A^ = 10^; Mean field 
case) . 



FIG. 7. (Color online) Relaxation time r as a funtion of 
8. Clearly r diverges as 6 ^ 6^. (Inset) Numerical fit to 
T ^\e- 6'c^^ with z ~ 1.9. (iV = 10*, ID case). 
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FIG. 6. (Color online) Variation of O, the average number 
of particles below the threshold energy 6 in the steady state 
{t > t), against threshold energy 9, following dynamics of 
Eq. P for ID. (iV = lO"). (Inset) Shows scaling fit (61 - Sc)" 
with /3 ~ 0.41. 



FIG. 8. (Color online) Variation of 0(t) versus t for different 
9 values for ID case (iV = 10*). 



the effective temperature is varied, has been studied [17| . 
Though the behavior is similar, the effective critical be- 
havior (exponent values) seem to be quite different. 

D. Results: Two dimensional model 

For the 2D lattice version, the particles are arranged 
on a square lattice, and this time one of the four near- 
est neighbors of i is chosen randomly as particle j. If 
the energy of the particle is below a prescribed thresh- 
old energy 9, then it collides with any one randomly 



chosen nearest neighbor j which can have any energy 
whatsoever, and the two particles will exchange energy 
according to Eq. ([!]). After each such successful colli- 
sion, the time is incremented by unity. This process 
is continued until steady state is reached. Variation of 
the steady state order parameter O against theshold 9 
is shown in Fig. [9l with exponent /? ~ 0.67 ± 0.01 and 
9c ~ 0.675 ± 0.005. Also, we find z ~ 1.2 ± 0.01 (see 
Fig. [TU]) and (5 ~ 0.43 ± 0.02 (see Fig. [Jill- 
All these estimated values of the critical exponents /3, 
and 5 are summarized in Table HI 
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This Model 


Manna Model 




ID 
2D 
MF 


0.41 ± 0.02 
0.67 ± 0.01 
0.97 ± 0.01 


0.382 ± 0.019 
0.639 ± 0.009 
1 


z 


ID 
2D 
MF 


1.9 ± 0.05 
1.2 ± 0.01 
0.83 ± 0.01 


1.876 ± 0.135 
1.22 ± 0.029 
1 


5 


ID 
2D 
MF 


0.19 ± 0.01 
0.43 ± 0.02 
0.93 ± 0.01 


0.141 ± 0.024 
0.419 ± 0.015 
1 



TABLE I. Comparison of critical exponents of this model with 
those of the Manna model [lal. 




FIG. 9. (Color online) Variation of O, the average number 
of particles below the threshold energy 6 in the steady state 
(i > r), against threshold energy 6, following dynamics of 
Eq. IT]) for 2D case. The simulation is done for lattice size 
1000 X 1000. (Inset) Shows scaling fit [O-e^f with /3 ~ 0.67. 
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FIG. 10. (Color online) Relaxation time (r) diverges as 6 
approches 6^ from both sides. (Inset) Numerical fit to r ~ 
\e - ec\~' for e < Oc { z ~ 1.2 ± 0.01). (simulations for 100 
X 100 system, 2D case). 
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III. FINITE SIZE EFFECT 

We have also studied the time variation of O at dif- 
ferent sizes {N) at 9c. Plots of 0{t)t^ as a function of 
t/TV^ for different values of the system size N (at critical 
point) are expected to collapse on a single curve. How- 
ever, as we have used a special dynamics vifhich never 
allo"ws the system to fall in the absorbing state, in this 
case the activity saturates at a small steady value (see 
Fig. [T2)) instead of shoviring the finite size cut-off. 

To study finite size effects in the decay of activity at 
the critical point, one has to remove the perturbation and 
allo'w the system to be trapped in the absorbing states. 
Using this dynamics, we have studied the effects of finite 
system size. Fig. [T^l show the decay of 0{t) with t at the 
critical point for MF, ID, and 2D systems, respectively. 
The insets in the corresponding figures show the data 
collapse. The fitting values of the exponent a are a = 
0.53 ± 0.02, 1.53 ± 0.05 and 1.55 ± 0.05 for MF, ID, and 
2D respectively, with 5 values given in Table-I. 



FIG. 11. (Color online) "Variation of 0(t) versus time t, shown 
for different 6 values for 2D case. At critical value 9c, order 
parameter follows a power law decay with exponent 5 ~ 0.43. 
The simulation is done for lattice size 500 x 500. 



IV. SUMMARY & DISCUSSION 

Inspired by the success of the kinetic exchange models 
of market dynamics (see e.g., 0-11]) and the observation 
that the poor or economically backwards in the society 
take major initiative in the market dynamics (see e.g., 
[H [3 and also [H [H), we consider an ideal-gas-like 
model of gas (or market) where at least one of the parti- 
cles (or agents) has energy (money) x less than a thresh- 
old (poverty line) value 9 takes the initiative to scatter 
(trade) with any other particle (agent) in the system, 
following energy (money) conserving random processes 
(following Eq. ([T])). For 6* — oo, the model reduces to 
the kinetic model of ideal gas with Gibbs distribution. 
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FIG. 12. Study of finite size effect: The variation of 0{t), versus t at 9 — 6c for systems of different sizes A'^ are shown. The 
plots with bare symbols correspond to the case when there is no absorbing state and those with symbols connected by solid 
lines correspond to the presence of absorbing state. The Figs, (a), (b) and (c) correspond to MF, ID and 2D respectively. 
Inset: The plot of Ot^ vs. t/N"' in presence of absorbing state, for different system size A'^ collapse onto a single curve for 
S = 0.93,0.19, and 0.43 and a = 0.53, 1.53, and 1.55 for MF, ID, and 2D respectively. 



The steady state is found to be ergodic (steady state re- 
sults are independent of initial conditions) . The pertur- 
bation employed in the frozen cases {xi > 9 for all i) also 
does not affect significantly the thermodynamic quanti- 
ties (e.g., the steady state value of O for 9 < 9c goes to 
with 1/N, as can be seen from inset (a) of Fig. [5]). 

In general, we find that the steady state distribution 
p{x) (see Fig. [1] for the mean field, where each parti- 
cle can interact irrespective of their distance from the 
active particle or agent having energy or money less 
than threshold or poverty line 9) differs in form signif- 
icantly from the Gibb's distribution, for finite values of 
9. The order parameter O, giving the average fraction 
of particles (or agents) having energy (or money) be- 
low 9 shows a phase transition behavior: = for 
9 < 9c and O 7^ for 6* > 9c (see Fig. [2] for mean 
field. Fig. HI for ID, and Fig. [9] for 2D, respectively). 
The critical values are given by 9c — 0.61,0.81,0.68 for 
mean field, ID and 2D cases, respectively. The varia- 
tion of O near 9c is quite universal (see Fig. [5]). We find 
O - (61 - 9c)'^ with /3 ~ 0.97, 0.41, 0.67 for mean field, ID 
and 2D cases, respectively. We also find that the relax- 
ation time T diverges strongly near 9c as t ^ {9 — 9c)~^ 
with z ~ 0.83, 1.9, 1.2 for mean field, ID and 2D cases, 
respectively. Finally, at 6* = 9c, 0{t) ~ where 
5 ~ 0.93,0.19,0.43 for mean field, ID and 2D cases, re- 
spectively (see Table H]). 

It might be mentioned here that the above exponent 
values are indeed very close to those of the Manna Uni- 
versality (MU) class (dl]; see also [2l|,[22]) in mean field, 
ID, 2D cases: Our estimates for (3 ~ 0.97, 0.41 and 0.67 
for mean field, ID and 2D cases, respectively, are quite 
close to /3 ~ 1, 0.38 and 0.64 for corresponding MU cases; 
5 ~ 0.93, 0.19 and 0.43 for mean field, ID and 2D cases. 



are also close to (5 ~ 1, 0.14 and 0.42 in the correspond- 
ing MU cases. However, it may be noted that signifi- 
cant differences in the above estimates do exist. Also, 
z ~ 0.83, 1.9, 1.2 for mean field, ID and 2D cases might 
be compared with z ~ 1, 1.87, 1.22 in the corresponding 
MU cases. These discrepancies could be due to finite size 
effect, and in that case the critical behavior of our model 
would belong to the MU class. As one can see, the esti- 
mated values of the exponents /?, 5 and z fit reasonably 
with the scaling relation 5 = P/z within our limits of 
accuracy. In this connection, it is worth mentioning that 
the violation of the above scaling relation has also been 
observed f23|, though such discrepancies seem to get re- 
moved if one uses all sample averages instead of averages 
over surviving samples [24| . In our case, however, this 
scaling relation seems to hold, as our simulation results 
correspond to all sample averages. 

In summary, when the energy threshold 9 is introduced 
in the kinetic theory of ideal gas such that the stochastic 
energy conserving scatterings between any two particles 
can take place only when one has energy less than 0, the 
gas system shows an intriguing dynamic phase transition 
at 9 — 9c, having the exponent values in the mean field 
(long range scattering exchange) , one dimension and two 
dimension, as estimated here using Monte Carlo simula- 
tion, are given in Table I. 
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